Equation-of-Motion Approach to Dynamical Mean Field Theory 
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We propose using an equation-of-motion approach as an impurity solver for dynamical mean field 
theory. As an illustration of this technique, we consider a fmite-C/ Hubbard model defined on the 
Bethe lattice with infinite connectivity at arbitrary filling. Depending on the filling, the spectra that 
is obtained exhibits a quasiparticle peak, and lower and upper Hubbard bands as typical features 
of strongly correlated materials. The results are also compared and in good agreement with exact 
diagonalization. We also find a different picture of the spectral weight transfer than the iterative 
perturbation theory. 
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Theoretical understanding of strongly correlated elec- 
tron systems has been a central and challenging problem 
in condensed matterphysics. Dynamical mean field the- 
ory (DMFT) QllQ has proven to be a very powerful 
theoretical tool that is not only tractable but also flexible 
enough to include material-specific details into the calcu- 
lations. The essence of the DMFT is to map the lattice 
problem onto a single impurity problem plus a set of self- 
consistency conditions through which the characteristics 
of a specific lattice model enters. In particular, all mod- 
els with only on-site Hubbard interaction can be mapped 
onto the single impurity Anderson model (SI AM) 0. 
The problem is then to find a suitable solver for the 
SIAM. For the last few years, several numerical and ana- 
lytical techniques to the SIAM have been explored. They 
include quantum Monte Carlo ( QMC) 0, EI j iterative 
perturbation theory (IPT) [1, 0, El El j an d exact diago- 
nalization (ED) |l2|,|l3j. Each method has its limitations. 
Except for some special models (l^j . QMC has difficulty 
in accessing very low temperatures, where statistical and 
finite time-step discretization errors become significant. 
Exact diagonalization, although exact, as the name im- 
plies, is limited by available computer time to a small 
number of orbitals, and hence it is difficult to obtain a 
smooth density of states. Iterative perturbation theory 
overcomes the drawbacks of these other two methods. 
However, in its original form, IPT is only valid at half 
filling. Away from half filling, it involves an ansatz and 
has instability problems. 

The aim of this paper is to present an equation-of- 
motion (EOM) approach that can overcome the difficul- 
ties of QMC and ED while still maintaining the numerical 
efficiency and other advantages of IPT. 

For concreteness, we demonstrate the method on the 
single-band Hubbard model for a Bethe lattice of infinite 
connectivity, z — > oo. The Hamiltonian is 



H 



(1) 



u, fl- 



and U is the on-site Hubbard interaction. 

We study this model using the DMFT approach, which 
maps the model onto a self-consistent SIAM as given by 
the Hamiltonian: 

#SIAM = YL tkcrC k<T Ck ° +^2( V k<? c k* d ° + V ka d Uka) 

+ ^2 tdcrd\d a + U n d ^n dl , (2) 



where ct (cfe CT ) are the creation (annihilation) operators 
for the conduction spin-i fermionic bath with the energy 
dispersion and d\ (d G ) for the impurity orbital with 
energy tda, U is the on-site Coulomb interaction between 
the electrons on the impurity, and Vka is the coupling 
between the bath and impurity orbital. The effective 
parameters and V^a enter the hybridization function 



A a (lUJ n ) = } 



(3) 



which is found by solving a self-consistent condition. 
Here u) n = (2n + l)irT is the Matsubara frequency 
for fermions with n being integer and T the tempera- 
ture. The impurity Green's function maps onto the site- 
diagonal Green's function of the original lattice Hamilto- 
nian by the relationship: 



SB) 



(e)de 



(e - n) 



where t (fixed) is the hopping integral for electrons be- 
tween z neighbors, rij the occupation number on site i, 



(4) 

where Gdaii^n) is the Fourier transform of the impurity 
Green's function G ti<T (r - r') = -(T T [d CT (r)4(r')]) de- 
termined from the impurity Hamiltonian Eq. J5J while 
Guo-(iu) n ) is the site-diagonal Green's function corre- 
sponding to the lattice Hamiltonian in Eq. QJ, pi B \e) is 
the Bloch density of states (the {7 = density of states of 
Eq. (1)) and \x is the chemical potential. This procedure 
indicates that A CT in the present context, in contrast to 
the original SIAM model [4j, is not known a priori. For 
infinite dimensional lattices or the Bethe lattice with in- 
finite coordination number as we are considering here, 
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the self-energy is site-diagonal and the DMFT becomes 
exact. In addition, for z — > oo, the non- interacting Bloch 
density of states is semiclliptic, which reduces the self- 
consistency condition to: 

A a (iu n ) = t 2 G da {iuj n ) , (5) 

with eda = — fJ- for the paramagnetic case. Keeping 
the hybridization function spin dependent allows the ap- 
proach to also be applied to the ferromagnetic ordering 
case. 

If one is able to solve the SIAM for arbitrary parame- 
ters (i.e., tker and Vka via A^), the lattice model can be 
solved through self-consistent iteration, by using the fol- 
lowing procedure: First, given a hybridization function, 
one solves the SIAM for the impurity Green's function (or 
self-energy); Second, through the self-consistency condi- 
tion, Eq. ||SJ), a new hybridization function is obtained. 
The procedure is then iterated until the self-consistency 
is achieved. 

The EOM approach corresponds to a resummation 
of low-order hopping process, but requires a decoupling 



scheme. We have used the one introduced by Appclbaum 
and Penn |l5j and Lacroix 

(APL), which is known 
to capture the right qualitative feature of physics at low 
temperatures. Briefly, the EOM method consists of dif- 
ferentiating the Green's function Gd^r — t') with respect 
to the imaginary time r, thereby generating higher-order 
Green's functions which must be approximated in order 
to close the equations for Gda(T — t'). Since the time 
derivative of a Heisenberg operator is determined by its 
commutator with the Hamiltonian, it follows: 

(-^-e drT JG da (T-r') = S(T-T') + UG drfa (T-T') 

+ J2v h : a G ka (T - t>) . (6) 

k 

On the right-hand side of Eq. Jfjjl, there are two new 
Green's functions: G ka {r - r') = -(T r [c fcCT (r)4(r')]), 
and the two-particle Green's function G dSa = 
-(T T [n dff (r)d CT (T)4 (>')])■ % usin S thc APL decoupling 
scheme to truncate this EOM one obtains: 



Gda(iu n ) 



1 - (n ds ) 



(do 



£/£, 



-e da -U-A a 



(rider) 



(da 



U 



us 2 



(7) 



where the occupancy probability is 

(nda) = J def FD (e)p IJ (e) , (8) 
with the spectral density p a {() = — lm[G da (e + ir])]/^ 



k 



where A™ = f FD (e k& ), A l g = 1 - /fdM, A<® = 1 
with the Fermi distribution function /fd(c) = [1 + 
exp(e/T)] _1 . By performing an analytic continuation, 
i.e., replacing iuj n with tu + irj, one arrives at the retarded 
Green's function, which is exactly identical to that de- 
rived directly from the real time approach 0, 0, 0] . 
The impurity Green's function Gda, as given by Eq. J7J), 
is exact both in the non-interacting limit (U = 0) and 
in the atomic limit (Vk a = 0). Since the accuracy of 
thc EOM method depends on how good the decoupling 
scheme is, this is an important confirmation on the qual- 



(tj is an infinitesimal). The hybridization function A CT 
represents the hopping of the electron with spin projec- 
tion a onto or out of the impurity site, while the other 
self-energies are due to the hopping of the electron with 
opposite spin projection a, 



,* = 1,2,3, (9) 
I 

ity of the APL decoupling. It is also useful to note that at 
low temperatures, the EOM solution based on the APL 
decoupling scheme has the expected Kondo resonance at 
the Fermi energy [Isj . originating from the self-energy, 
Sict, which is due to the virtual intermediate process 
where the impurity orbital is occupied by an electron 
of opposite spin a. 

Physical quantities like the spectral function require a 
knowledge of thc Green's function on the real frequency 
axis. For IPT, direct calculations of thc self-energy on 
the real frequency axis are very complicated because of 



^da + tdo — e kS 



1 

(da — (da — U + 6ka 



2 
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multiple-dimension integrals. Instead one first performs 
the calculation on the imaginary frequency axis and then 
computes the spectral function by performing an analytic 
continuation by a Pade approximation. For this kind of 
treatment to be successful the quantities involved in the 
Green's functions and the self-energies need to be ex- 
pressed in terms of the hybridization function as a whole. 
In the current approach, the three self-energies in Eq. 



J 



cannot be represented as functions of A CT . However, by 
introducing a line-width function 



r CT (cj) = Im[A CT (w + iij)] , (10) 

7T 



one can instead rewrite Eq. as: 



Hia(iw n ) = / deAW(e)r ff (c) 



tds — e —iu>„ 



U-6 



1,2,3 



(11) 
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FIG. 1: The spectral density as a function of energy for [7 = 4 
and the hole doping 8h — (a) and 0.12 (b). Thick solid line: 
EOM approach. Thin dashed line: ED (6 sites). For clarity, 
the intensity from the ED is scaled down by a factor of 16 (a) 
and 8 (b). The parameter 77 = 0.01. u) — corresponds to the 
Fermi energy. 

In order to work on the imaginary frequency axis, at 
each iteration step a Pade approximation for T a would 
be required, which would introduce significant numerical 
errors and make the convergence difficult. Since the occu- 
pancy probability given by Eq. © and the self-energies 
given by Eq. I|ll|) all involve only one-dimensional in- 
tegrals, we find it is more efficient to implement the 
technique directly on the real-frequency axis. For this 
purpose, we divide the relevant energy range into seg- 
ments and use the Simpson's rule in combination with 



the Neville's polynomial interpolation algorithm to do 
these integrals [l!|. In the actual calculation, we use 
two iteration loops: The inner loop is for the occupancy 
propability (rida) for a fixed line- width function T a . Af- 
ter the convergence of the inner loop is achieved, T a is 
updated through Eqs. JSJ) and i|10fl in the outer loop. 
We find the algorithm is very efficient: in most cases a 
solution is found within 20 iterations. It is also remark- 
able that the algorithm provides stable convergence with 
respect to different guesses for the initial hybridization 
function used to start the calculation. 

In the following discussions, all energies are measured 
in units of the half bandwidth D = 2t = 1 . We take the 
temperature T = 0.01. In Fig. we show the spectral 
density, p a (u), as a function of energy for U — 4 and for 
values of the hole doping 8u (= 1 — S<j( ndcr )) °1 an< l 
0.12. As shown in Fig. Q^a), in the undoped case (i.e., 
half filling), the spectral density displays only the lower 
and upper Hubbard bands separated by a Mott insulator 
gap. Due to the particle-hole symmetry, these two bands 
are symmetric to each other around the Fermi energy. 
At finite doping (e.g., 6 = 0.12 as shown in Fig. QJb)), 
the spectral density exhibits a sharp quasiparticle reso- 
nance peak at the Fermi energy, and two broad satellite 
features corresponding to the lower and upper Hubbard 
bands, as would be expected for a doped Mott insulator. 
Since the particle-hole symmetry is broken upon doping, 
the spectral density is no longer symmetric with respect 
to the Fermi energy. To demonstrate the accuracy of the 
present approach, we have also run the ED algorithm 
developed by Caffarel and Krauth 0] f° r 6 sites with 
the same parameter values. The output is shown by thin 
dashed lines in Fig. ^ As one can see, the results ob- 
tained from the EOM approach compare well with those 
from the exact diagonalization method. Like IPT, in ad- 
dition to computational efficiency, the advantage of the 
EOM approach is that it can give a very smooth shape 
for the spectral density. The ED method can at best give 
a correct overall spectral density distribution with sharp 
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spike-like structures due to the small number of orbitals 
being used. Because the use of a large number of orbtals 
in the ED method would exceed current computational 
capabilities, we conclude that the EOM approach may 
be a useful alternative. 




FIG. 2: (Color) The evolution of the spectral density with 
increasing hole doping S h = 0.06, 0.16, 0.26, 0.44, 0.64, 0.88 
(counted from the left). The parameter r\ = 0.02. 

With increased confidence due to our success in calcu- 
lating the spectral function, we now turn to the evolution 
of the spectral density with hole doping. These results 
are shown in Fig. [21 The most striking features are as 
follows: When the chemical potential crosses the edge of 
the lower Hubbard band (i.e., when hole doping begins), 
a Kondo resonance peak shows up at the Fermi energy. 
At the same time, an anti-resonance dip occurs in the up- 
per band. As the doping increasing, there continues to be 
a peak-like feature near the Fermi energy that gradually 
decreases in intensity. The effects on the upper Hubbard 
band are even more dramatic, as the anti-resonance dip 
evolves into a feature that shifts more deeply into the 
upper Hubbard band and eventually becomes more of a 
plateau, with the spectral weight continuing to be trans- 
ferred out of the upper Hubbard band. The hole doping 
dependence of the spectral weight reflects the different 
parameter regimes of the SIAM. Within the EOM ap- 
proach, one does not find a narrowing of the upper Hub- 
bard band. This picture of the spectral density evolution 
is different from that based on the iterative perturba- 
tion theory. Since the two approaches are so different, 
we find it very difficult to draw a judgement on which 
picture is more reasonable without a large-size ED cal- 
culation to check the results. We have also studied the 
electron doping case, and found a similar spectral weight 
transfer (now out of the lower Hubbard band). 

To summarize, we have proposed an EOM approach as 
an impurity solver for dynamical mean field theory. The 
solution for the impurity Green's function is exact in the 
non-interacting and isolated-site limits. The technique 
has been implemented directly on the real- frequency axis, 
which turns out to be computationally efficient. As an 
illustration, we have studied a finite-?/ Hubbard model 



defined on the Bethe lattice with infinite connectivity 
at arbitrary filling. Depending on the filling, all typi- 
cal features of strongly correlated materials such as the 
quasiparticle peak at the Fermi energy, and lower and 
upper Hubbard bands, have been produced for the spec- 
tral density. The main results have also been compared 
with exact diagonalization and found in good agreement. 
In addition, a new picture of the spectral density evolu- 
tion as a function of filling factor is found, which should 
be tested by future exact diagonalization studies with a 
large number of orbitals. The extension of this approach 
to the multiorbital case constitutes a future investigation, 
which will be very useful for ab initio DMFT applications 
to real materials. 

Note Added: In the final stage of this work, we 
became aware of a preprint by H. O. Jeschke and G. 
Kotliar |20j |. where a similar approach was proposed. 
Their study focused on the case of an infinite Hubbard 
interaction such that the technique could be directly im- 
plemented in Matsubra frequency. In that case, only the 
quasiparticle peak and lower Hubbard band are the rel- 
evant characteristics. It is encouraging to see that in 
the comparable region, both of our studies show a very 
similar shape for the spectral density. 
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